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Summary 



In order to choose a numerical method for solving the time dependent equations of radiative transport, we obtain 
an exact solution for the time dependent radiation field in a one dimensional infinite medium with monochromatic, 
isotropic scattering for sources with an arbitrary spatial distribution and an arbitrary time variation of their power. 
The Lax-Wendroff method seems to be the most suitable. Because it is assumed that radiation delay is caused by 
the finite speed of light, the following difficulty arises when the numerical method is used: the region of variation of 
the variables (dimcnsionless coordinate r and time t) is triangular (the inequality r < t). This difficulty is overcome 
by expanding the unknown functions in series in terms of small values of the time and coordinate. By comparing the 
numerical and exact solutions for a point source with a given time dependence for its power and with pure scattering, 
the steps in the variables required to obtain a desired accuracy are estimated. This numerical method can be used to 
calculate the intensity and polarization of the radiation from sources in the early universe during epochs close to the 
recombination epoch. 

Keywords: radiation scattering: time dependent 

1. Introduction. Analytical and numerical solutions are given in this paper for a simple problem of time dependent, 
monochromatic scattering of radiation in a so-called one dimensional medium. 

An idealized medium is refered to as one dimensional if a photon, either preserves the direction of its motion upon 
scattering in this medium, or its direction of motion reverses (backward-forward scattering), so that it propagates 
along a straight line. The assumption of a one dimensional medium is equivalent to assuming that the scattering takes 
place with an indicatrix equal to the sum of spike (delta- function) indicatrices. This kind of scattering problems was 
examined in the early development of radiative transport theory [1, 2, 3, 4, 5, 6]. The description of scattering in an 
ordinary (three dimensional) medium reduces approximately to the transport equations for a one dimensional medium 



Here we examine the propagation of radiation in a one dimensional medium with invariant optical properties. It 
is assumed that the scattering is isotropic, monochromatic and, as usually assumed for monochromatic scattering, 
such that the time delay of the photons is caused by the time they spend in transit, while scattering events are 
instantaneous. The sources of primary radiation can be nonisotropic, with an arbitrary spatial distribution, and have 
an arbitrary time dependence. 

Since the medium is stationary, we can introduce the optical path t, measured from some point in units of the 
mean free path of the photons. We restrict ourselves to scattering in an infinite medium, where the optical path r 
ranges from — oo to oo. The time t is also measured in units of the mean time between collisions. Sources can act from 
time t — — oo , when there is no radiation field. The solution of this problem is obtained here in an explicit analytical 
form. 

The main object of attention is a point source with a given time dependence of its radiated power. The existence of 
an exact solution has made it possible to choose a numerical method for solving equations similar to those for scattering 
in a one dimensional medium. It is planned to use this numerical method in a separate article for calculating the 
evolution of the intensity and polarization of radiation from sources in epochs of the early universe close to the 
recombination epoch (see [7] and [8], for example). 

2. Basic equations. We denote the intensity of the radiation propagating in the direction of increasing and 
decreasing optical depths by r) and 7_(£,t), respectively. Then the two equations for the evolution of the 
radiation field can be written in the form 



[!]• 



j±(t,r)±4(t,r) 



I±{t,T)+B±(t,r). 



(1) 
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Here a time derivative is indicated by an overhead dot and a derivative with respect to the optical path length, by a 
prime. The two source functions, as always, consist of two parts, one characterizing the power of the primary source 
and the other describing the scattering: 



B±(t,r) =g±(t,T) + XJ c (t,r), 



(2) 



where A is the probability of survival of a photon in each scattering event. Here we have introduced the average 
intensity J c {t 7 T). We can immediately introduce the radiative flux, as well: 



I+(t,T)+I_(t,T) I+(t,T) - I-{t,T) 
Jc(t,T) = , H c (t,T) = 



(3) 



2 ' 2 

Adding and subtracting Eqs. (1) yield equation for the average intensity and the flux: 

j c (t,r) + n' c (t,r,t) + (i - x)j c (r,t) = /j(T,t), n c (t,r) + j^(t,r) + n c (t,r) = / H (t,r), (4) 

where 

n , (t. n-\ 4- n (+ -A n , (t. t\ - n (i : t\ 

(5) 



m,r) = 9+{t ' T)+ J- {t ' T) , Mt,T) = 9+{t > T) - g - {t ' T) . 



3. Laplace transform. An efficient way of solving time dependent problems is to use Laplace transforms with 
respect to time. We shall indicate the transform with a tilde over the transformed function: 



oo oo 

J(t,8)= J J c {t,T)e~ st dt, H{t,s)= J H c {t,r) 



e~ st dt. 



Applying the transform to Eqs. (4) for zero initial conditions (at t = — oo) yields 

H\t, s) = -(s + 1 - X)J(t, s) + /j(r, s), J'(t, s) = -(s + l)H(r, s) + / H (r, s). 



(6) 



(7) 



Let us find a solution to Eqs. (7) assuming that the transform parameter s is real and non-negative. Finding 
the general solution of the homogeneous equation and a particular solution for the non-homogeneous equations, while 
noting that the solutions must be finite at r — > ±oo, we obtain 



oo 



H(t,s) = 



-fe|r-n| 



S + 1 ~ 

—, — /j(ti, s) + /h(ti, s) sgn(r - n) 



/j(n, s) sgn(T - n) H — / H (ri, s) 

s + I 



dn, 



dri. 



(8) 
(9) 



Here the positive root k — a/(s + l)(s + 1 — A). 

4. Inversion of the transforms. In order to invert the Laplace transforms, we use Eqs. 4.17 (5)-(9) of the Bateman 
and Erdelyi handbook [9] . Excluding the singular solutions which arise on inverting the parts of the transforms which 
behave as e~( s+1 ~ A / 2 )l T ~ Tl l, we write the result of the inversion in the form 

oo 

J c (t,r) = \ J dne-^-^lr-Til [/j(t- | r -n|,Ti) + M* - k-n|,ri)sgn(T-Ti)] + 

— OO 

oo t-\r-T\\ 

+ \ J dn J di 1 e- (1 - A)(t -* l) [Gj.T(A(t-t 1 ),A(T-T 1 ))/ J (t 1 ,T 1 )+G.TH(A(t-t 1 ),A(r-T 1 ))/ H (ti,r 1 )],(10) 

— oo —oo 

oo 

n c (t,r) = \ j dne-WK^ [fj(t-\r-r 1 \,T 1 )sgn(T-T 1 ) + f K (t-\r-r 1 \,T 1 )} + 

— oo 

OO t—\T—Tl\ 

+ I / ^ / d *i e ~ (1 " A)( '^ l) [ GHJ ( A ( < - < i)' A ( T - T i))/j( < i' Tl )+ GHH ( A ( t -*i)' A ( T - T i))/ H (*i' Tl ^ 
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This solution is expressed in terms of the Green functions 



G.T.T(t,r) = e-*/ 2 
G H H(*,r)= e -*/ 2 



w 
t 

IV 



G m (t,r) = Ghj(t,*) = e-*/ 2 -/! (^) . 

w \ 2 J 



(12) 

(13) 
(14) 



In these formulas I n (z) is the Bessel function of imaginary argument and w = \lt 2 — r 2 . The Green functions have 
meaning for t > |r|, which sets the upper limit of integration on t x in Eqs. (10)-(11). The values of the functions at 
the boundary of the region in which they are defined, i.e., t = |r|, are 



Gjj(|r|,T)=e- 



M/2 



G 



HH 



|r|,r)=e-M/ a (ili-l) 1 G aH (|r| ) r)=e-M/ 2 ^. 



It can be verified by direct substitution that the functions (12)-(14) satisfy the homogeneous equations 
Gjj + G' H j = 0, Gjh + G^jh = 0) Ghh + Ghh + G' H} — 0, Gjh + Gjh + G JT = 0, 



(15) 



(16) 



which can be used to confirm that Eqs. (10)— (11) do, indeed, determine the solutions of Eqs. (4). 
Function (12) was obtained previously by Minin [4]. 

5. Point source. Let us consider the important special case of a time dependent, isotropic source acting at t > 0. 
The following functions apply in this case: 



/j(t,r) =S(r)C(t), / H (t,r)=0, 



(17) 



where the function £(t) characterizes the variation of the source power with time. We shall assume that it is continuous 
and falls off rapidly enough so that the effective duration t s of the source is finite. The corresponding average intensity 
and flux are single integrals. For t> |r|, 



Jc(t,r) - ±e-V-VQ\T\£{t - |r|) + J ; (i,r), 
H c (t, t) = hr^-WWCit - \t\) sgn(r) + H^t, r), 



(18) 
(19) 



where the integral terms are given by 



(- — | T | t — T 

Ji(t,r) = ± J dt ie -( 1 - A )(*- tl )G JJ (A(i-i 1 ),Ar)£(t 1 ), Hi(t,T) = j J dhe-^-^-^G^iMt - h), Xr)C(h). 
o o 

(20) 

The average intensity is continuous for all r and t > |r|. For t = and C(t) > 0, i.e., at the site of the point source 
and during the time it acts, the flux has a discontinuity: H(+0, t) — H(— 0, t) = C(t). At the points t = ±t, where the 
medium only begins to radiate, the solutions have the following limiting values: 



J(\t\,' 



1 -(i-A/2)|r| r(0)j H (\T\,r)=J(\T\,r) S gn(r). 



(21) 



It is easy to determine how the solutions behave at late times compared to the time the source is active and at 
large distances from it, i.e., for t » t s and f>r. We shall assume that the ratios t s /t and r/t are small, while the 
ratio to = r 2 /t is of the order of unity (i.e., t /t is small). For A = 1, we have 



Ji(i,r) 



-to/4 



2Vnt 



i ii 



to 



1 



16 

4 



2 4 

f3 



(tl 11 2 21 3\ / 4 4 27 2 3 3 \1 1 1 
+Cl 1 64 _ 32*° + 16*° " 8 ) + C ° [512 16 + 64*° " 8 <0 " 32 ) \ T 2 ) 



(22) 
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In these equations the moments of the source power are 

oo 

C n = J C(t)t n dt, (24) 
o 

which are assumed to be finite. 

The main terms in the asymptotes are solutions of the diffusion equation, and a diffusive relation holds between 
the average intensity and flux, i.e., H(t,r) = — J'(t, r). The flux falls off with time substantially more rapidly than 
the average intensity. 

6. Direct and scattered radiation. The solutions of Eq. (4) for A = represent the radiation propagating directly 
from the source, without undergoing any scattering. For a point source, these solutions are 

J.(t,r) = ie-l r l£(t- |r|), fr„(t,T) = J„(t,T)sgn(T). (25) 

The scattered (diffuse) radiation is characterized by quantities that go to zero for A = and are equal to the differences 
of the corresponding functions (f8)-(19) and (25), i.e., 

J d {t,r) = J s {t,r) + Ji(t,r), H d {t,T) = J s (t,T)sgn( T ) + H i (T,t), (26) 
where the expression outside the integrals is 

J s (t, t) = A(r)C(t - |r|), A(r) = \ (e^-V^M _ e -|r|) . (27) 

The functions J d {t,r) and H d (t,r) obey the equations 

j d {t, t) + H' d (t, t) + (1 - X)J d (t, t) = J.(t, r), H d (t, t) + J' A (t, t) + H d (t, t) = 0. (28) 

The boundary conditions for these functions are analogous to Eq. (21), but instead a single exponent it is necessary 
to take the difference, i.e., to replace the factor in £(0) by the function A{t). At the site of the source (for r = 0) the 
average intensity and flux of the scattered radiation are continuous, with 

J s (t,0)=H d (t,0)=0. (29) 

The asymptotes of the functions Jd(t,r) and Hd(t,r) coincide with the asymptotes of the functions Ji(t,r) and 
Hi(t, t), since the expressions outside the integral fall off more rapidly than the integrals. 

Because of the obvious symmetry of the average intensity and the asymmetry of the flux with respect to r, in the 
following we consider r > 0, i.e., we consider the radiation field from one side of the source. 

7. Power-law sources. We now consider a rather general kind of source whose power can be expanded in a power 
series in time near the onset of their activity: 

oo 

C(t)=t»L(t), L{t) = Y,Lit\ (30) 

1=0 

where /i > 0, while Lq ^ 0. For source of this type it is appropriate to introduce the new unknown functions (with /i 
and A as parameters) 

Jd{t,T)=y*J(t,T), H d (t,T)=y»H(t,T), y = t-r. (31) 
These functions, again, obey the equations 

j + H' + ^(J-H) + (1- A)J= \e~ T L{y), H + J' + ^{H - J) + H = 0. (32) 



The boundary conditions for these functions are derived from Eq. (21): 



J(t,t) = H(t,t)=A(t)L . 
The integral terms near the boundary t — t for power-law sources behave as (y <C 1) 

y 



Ji(T + y,r) 



fl-i(r + y,T)~ | 



LqGhj (t,t) y 



G H j(t, t)Xi + 



ioG^(r) 



M + 1 /i + 2 y nJ v ' n+l 

Here the values of the Green functions at the boundary are given by Eq. (15), while 



(33) 

(34) 
(35) 

(36) 



8. Expansions for small t and t. The distinctive feature of the problem is that the region of variation of its 
parameters has the shape of an infinite triangle with boundaries in the (t, r) plane formed by the rays t = and r = t 
emerging from the coordinate origin. Near the vertex of this triangle it is extremely difficult to discretize the variables, 
so that it is impossible (or very difficult) to use a numerical method from the time of the very onset of the source, 
t = 0. Thus, it makes sense to move slightly away from the vertex point to some other time. Near this point the 
variables have small values and, if the function L(t) has several derivatives at t = 0, then all the functions, including 
the unknown functions, can be expanded in Taylor series in the neighborhood of this vertex. If L(t) is analytic, then 
expansions in power series are possible. Note that our functions arc not analytic with respect to the argument r, since 
they depend on |r|. Thus, the power law expansion, we shall obtain for t > 0, cannot be extended to negative r. 

It is easy to see that for a power law source whose power is given by Eqs. (30), because the equations are linear, it 
is possible to resolve components of the radiation field belonging to different powers in the expansion of the function 
L and proportional to the coefficients L t , i.e., one can write the average intensity and flux in the form 



J(t, r) = ^ y l U Jt(T, y), H(t, r)=J^ v'Wr, y), 



(37) 



1=0 



1=0 



where y — t — r. 

Substituting Eqs. (37) in Eqs. (32), we find equations for the functions which are the coefficients of the expansions: 



Ji + Hi + (1 - A) J; + ^(J ; - Hi) = Kr T , Hi + J[ + H- ^(Ji - Hi) = 0. 

y ^ y 



(38) 



We now transform to the variables r and y, on which these functions depend. Then the new equations take the form 



"97 ■+ d-y {Jl Hl > '+ y 



' ^(Jt-HO + a-XV^^e- 



^- ^{Ji Hi) -^{Ji- HO +Hi=0. (39) 
or ay y 



Adding these equations yields 



^-(Ji + H t ) + (1 - X)J t + Hi = 



(40) 



As y — > 0, Eqs. (39) imply that J;(r, 0) = Hi(t,0), while in the expansions (37) only the zeroth order terms 
Jo(t, 0) = Ho(t, 0) remain. Taking the same limit in Eq. (40) and using the equality above yields the following, easily 
integrable linear differential equation: 



dJo(r,0) 
dr 



+ 1 



^ J (r,0)^ e - 



(41) 



Solution of Eq. (41) gives the function that depends on r in the term under the integral in Eq. (27) subject to the 
boundary condition (33): Jo(t,0) = A(t). 
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It is clear from Eqs. (38) or (39) that the index /j and number I appear in them only as a sum. Thus it is 
appropriate to use the notation + 1 + 1 = z and regard this sum as an implicit argument of the functions Ji (r, y) and 
Hi (t, y) . Hence, we can write the expansion for these functions at early times and for small values of the coordinate 
in the form 



oo m 



(42) 



m— j — m— j— 

where the parameters A and z are indicated as arguments of the coefficients with their two indices. In the following 
we shall often omit these arguments. 

Equality of these functions at y = requires that J m ,o( z , A) — H m0 (z, A). The fact that the flux is zero at t = 
leads to the following condition 



Hi(0,y)= E#m, m (2,A)2/ m = 0, 



(43) 



m=0 



so that H m , m (z, A) = for all to. 

Substituting expansion (42) into Eq. (39), we obtain a pair of recurrence relations (omitting the arguments z and 

A), 



(m + 1 - j)H m+lij + (z + j)(J m+ ij + i - H m+ i ij+1 ) + (1 - A) J, 



A(-l)' 



(m + 1 - j)J m +i,j - (z + j)(J m +i,j+i ~ H m+ i ij+ i) 



H, 



m,3 



2 to! 
= 0. 



■$j,0, 



Adding these two equations yields a simple formula which follows from Eq. (41): 

A(-iy 



(m + 1 - j)(J m+1 j + H m+hi ) + (1 - X)J m ,j 



■Sj,o- 



The value j = gives 



from which we find 



(to + l)J m +i,o + 



1-^ 

2 



Jm : — ^ 



Jm,0 — H m fi — 



!(-!)' 



to! 



TO 



- 1 



(44) 
(45) 

(46) 
(47) 
(48) 



which is the coefficient in the expansion of A(t) . 
Setting to = and j = in Eqs. (44)-(45) gives 



#i,o + z(J 1A - H 1A ) + (1 - A) J ,o = 



A 



^i,o — z (^i,i — #1,1) + #o,o = 0. 



#1,0 — #1,1 



(49) 



A 



0, this leads to the value of yet another coefficient: Ji i = — . 

4z 



Given that Jo,o = #o,o = 0, Ji.o 

In general, for a given value of to, 2(to + 2) new coefficients J m +ij and H m+ ij appear, where j = 0, 1, to + 1, 
but only to 2(to + 1) equations. However, the values of J m +i,o = # m +i,o are known, while H m+ i, m+ i = 0, so that 
one of the equations is redundant. This equation, the sum of a pair of equations (for j = 0), i.e., Eq. (47), has been 
solved for all to. The procedure for solving the equations for all to is the same. First we find the difference in the 
coefficients with one value of j, and then their sum, going from smaller to larger j. Here z enters as a parameter. 

9. Expansions of the exact solutions. The method of obtaining recurrence relations for the coefficients of the 
expansion can be applied to more complicated and more general equations, including those with variable coefficients 
in front of derivatives. We suppose to use it for calculating the radiation fields in an expanding universe. In this work, 
it is possible to obtain expansions for the functions from their explicit expressions. 

Successively expanding the Bessel functions and exponents under the integral sign in Eqs. (20) and then integrating 
with respect to t\ and equating all the terms in the form of functions of r and y, for the average intensity we find the 
expansion coefficients for m > 1 and 1 < j < to to be 



[(m-i)/2\ 

. A ^r-^ — (TO — 2fc — 1)! 

Jm,j{z, A) = -^Bj(z) ^ Jm-l,kW ( m _2k-iV 
k=0 V ■"' 



(50) 
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where the square brackets in the upper limit of the sum denote taking the integer part, while 



j'-i 



B 3 (z) = Y— , ( 1 1 ) ,„ ; ., , Z = fi + l + l. 



(51) 



The coefficient under the summation sign, which depends on A, is given by 



Jm,fe(A) = - } ^ 



y2nyn — 2n 



n—k 



2 4n n\(n- k)\{m- 2n)\ \ 4Ai n+1 



A m — 2n 



(52) 



where Ai = 1 — A/2. 

The coefficients in the expansion for the flux are obtained in a similar way (m > 1, 1 < j < m): 



A 



[(m-j)/2] 



H m+1 ,j(z,X) = — Bj{z) Y H m-l,kW 



fe=0 



(m-2k- 1)1 
(m — 2k — j)\ ' 



-ffm,fc(A) — 



(— \) k+r ' 



i/2] 



}2n^m — 2n 



n—k 



2 in {n+ l)!(n - k)\(m- 2n)V 



(53) 
(54) 



The coefficients obtained from the exact formulas and from the recurrence relations (44)-(45) are the same. 



10. Bell-shaped sources. We now consider one concrete examples of a time dependence for the source power: 
sources acting for a finite time t s , delivering a total energy equal to 1, rising and turning off smoothly. For a source 
of this type, we can assume a time dependent profile of the form 



L(t) = c, 



hi 



t 

1 — cos I 2ir— 



L s 



t 



= 2c M ^- sin z [ 7T- ] , < t < t B , n > 0, 



(55) 



where c M is a normalizing coefficient. 

Table 1 . Estimated Accuracy of the Asymptotes at t = 20 



T 


J(t,r) 


Ja,s(t,T) 


H(t,T) 


H as (t,r) 





0.062695 


0.062687 


0.0000000 


0.0000000 


5 


0.046926 


0.046928 


0.0058817 


0.0058821 


10 


0.018610 


0.018589 


0.0049194 


0.0049089 


14 


0.048386 


0.048526 


0.0019490 


0.0019418 


18 


0.004726 


0.005010 


0.0002925 


0.0003137 



Because the source acts for a finite time, the upper limit in the integrals of Eqs. (20) should be taken to be 
min(t — |T|,t s ). The rest of the equations are valid without changes. In particular, the expansions for a bell-shaped 
source can be obtained by noting that in the expansion of Eq. (55), only even coefficients are nonzero: 



21 



L21 = 2L Q i ^— ( — ) , L21+1 = 0, 



(-1)' 



(2/ + 2)! \t s J ' 

The moments of the power are given by the series 

00 

c n = 471-% 



Lo = L(0) = 27r^. 



(27T) 



21 



(21 + 2)121 + 1 + fi + n' 



(56) 



(57) 



The normalizing coefficient is determined by the condition £ = 1. In particular, c = 0.11223, C\ = 0.41023, c 2 = 1 and 
c 3 = 2. 

As a comparison, Table 1 lists the asymptotic (Eqs. (22)-(23)) and numerically exact values of these quantities 
for t = 20 and a number of values of r. Naturally, the accuracy deteriorates with increasing r. 
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0.0 




2 4 6 8 10 i 

Fig. 1. The term outside the integrals in Eqs. (18)— (19) 
as a function of t for t s = 1, /U = 2, A = 1 and r = 0(1)9. 

Figure 1 contains 10 plots of the part of the average intensity and flux outside the integrals in Eqs. (26) (function 
(27)) for r > while figures 2 and 3 show the integrals (20) for A = 1, /j, = 2 and t s = 1 as functions of time for several 
values of r. The curves correspond to those values of r from which they begin on the abscissa. In the figures we can 
see only the tendency of approaching to the asymptotes. 

11. Numerical method. Several direct methods of solving the system of Eqs. (28) based on discretizing the 
equations were tried. The most suitable was the Lax-Wendroff method. This method is usually prescribed for partial 
differential equations with zero right hand sides [10]-[13], but it is easy to generalize it to systems with nonzero right 
hand sides. 

This method turned out to be sufficiently stable if it was applied to equations of the form (28) for a source with 
power given by Eq. (30) with integer values of the parameter /i. Here we illustrate this numerical method for the case 
of pure scattering, A = 1, and fi = 2. For brevity we omit the index d. The equations take the form 

j + H' = \e- T C{y), n + J' + 7i = 0. (58) 

If the unknown functions have already been found at some time t for a set of values of the coordinate r, then 
according to this method the values for subsequent times and the same coordinates are found using a second order 
Taylor expansion, 

J(t + At,T) = J(t, T ) + j(t,T)At + J(t,T)^-, H(t + At,r) =H(t,r) +H(t,r)At + H(t,T)^-. (59) 

The time derivatives can be expressed in terms of derivatives with respect to the coordinate using Eqs. (58). 

Let us discrctizc the variables, taking equal step h in time and in the coordinate: tj = hi, tj — hj, i = 0(l)io, 
j = 0(l)i. The region of the discretized points is shown schematically in Fig. 4. As the initial time we take 
t* = ti„ = i*h, where the number is defined below. We calculate the values of the unknown functions at this time 
using their expansions (37) and (42). The coefficients in the expansions are found either using the recurrence relations 
(44)-(45) or from their exact expressions (50)-(54). 

We denote the values of the unknown functions at the nodal points by J7i,j = J7(ti,i~j) and Hij = Ti{ti,Tj) and 
that of the source power, by £i-j = C(h(i — j)). The corresponding values of the derivatives are indicated, as before, 
by a dot and a prime. A prime also indicates a derivative with respect to the argument of the function C(t). The 
numerical values of the derivatives with respect to the coordinate are calculated using the formulas 

qi ^ —Ji,j-\ sjii ^ JiJ+l ~2>Ji,j + >Ji,j-l nil ^ + l njti ^ 

Jij ~ 2h ' c ~ //-' ' ~ 2h ' ' ~ h 2 ' { ' 
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2 4 6 8 10 12 14 t 
Fig. 2. The integral term in J[(t, r) of Eq. (20) as a function of t 
for t a = 1, \i = 2, A = 1 and r = 0(1)9. 
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Fig. 3. The integral term in Hi(t,r) of Eq. (20) as a function 
of i for i B = 1, p, = 2, A = 1 and r = 0(1)9. 
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^i* — fit* t% ti-\-i ^ 

Fig. 4. The scheme for discretizing the equations. 

As a result of discretization, we obtain a system of discrete recurrence relations of the type (for i > i t , with j = 
1, 1) 

J»U - \e~ ih [2£i-i + hC>_ 3 ] - (l - K < 1 H - > : + -W+^-i , (61) 
The boundary conditions (33) imply that 

Ji,i = Hi ti = (63) 

for all i. The values of the unknown functions below the boundary t — r (at the points with j = i — 1) indicated by 
asterisks in Fig. 4, were obtained by interpolation using Newton's formula with three points j = i + 1, i — 1, i — 2: 

Ji+IA = Ji+l,i-l — i7i+l,i-2/3, Wj+i^ = Hi+ij-i — Wj + i i i_2/3. (64) 

The flux at the source is equal to zero, and the value of J at r = was found by extrapolation. Thus, for all i, we 
have 

— = 0, JiS) = 3(j7i4 — ^,2) + <^i,3- (65) 

The following order of calculations was employed: first, for the chosen step size h the values of the unknown 
functions were calculated from their expansions for i — i*. Then a transition to larger values of i was made in 
succession: for j = 1, — 1 from i to i + 1, using Eqs. (61)-(62), for j = i + 1 using Eqs. (63), and for j = and 
j = i using the values already calculated according to Eqs. (64) and (65). 

During the calculations, with increasing i it was necessary to increase the step size, taking the calculated values 
for the last i as the initial values. The former were taken to be h = l/2 m , t* = l/8(i* = 2 m ~ 3 ), with m chosen to be 
m = 8 — 10, depending on t s . At t = 5 and then at t = 12, the step size was doubled. The calculation was continued 
up to some number i = i when the estimated functions attained their asymptotes with the same accuracy as the 
exact functions. 

A comparison of the numerically determined values of J(U,Tj) and H(ti, tj) with those calculated with the exact 
formulas showed that for m = 8 and t s = 1, the maximum relative error was 10~ 3 , while for t s = 2 it was 10~ 4 . 
The error increases as the parameter t s is reduced, since the function (55) becomes narrower and higher (closer to a 
(5-function). When the allowed accuracy is reached, it is necessary to reduce the step size. 

A similar numerical scheme, based on a predictor-correction method, has been developed by MacCormac [14]. It 
can also be used to solve these problems. 

12. Conclusion. The exact solution of the equations obtained here have made it possible to choose a method 
that is suitable for solving them numerically and can be used to estimate the parameters required to achieve a given 
accuracy. It is proposed that this method be used for calculating the evolution of the luminosity and polarization of 
sources in the universe during epochs close to the recombination epoch. 
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